# Partisan Collaboration in Policy Adoption: 
# An Experimental Study with Local Government Officials
# Yixin Liu 2024
# Policy Studies Journal
# Replication file for descriptive analysis

####Packages####
library(dplyr)
library(tidyverse)
library(ggplot2)
library(cowplot)
library(egg)
library(gridExtra)
library(grid)
library(ggpubr)
library(maps)
library(ggmap)
library(mapproj)
library(plotly)
library(pastecs)

####Functions####
# Share a legend between multiple plots
grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {
  
  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x + theme(legend.position="none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)
  
  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl),
                                            legend,
                                            ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend,
                                           ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))
  
  grid.newpage()
  grid.draw(combined)
  
  # return gtable invisibly
  invisible(combined)
  
}

element_textbox <- function(...) {
  el <- element_text(...)
  class(el) <- c("element_textbox", class(el))
  el
}

element_grob.element_textbox <- function(element, ...) {
  text_grob <- NextMethod()
  rect_grob <- element_grob(calc_element("strip.background", theme_bw()))
  
  ggplot2:::absoluteGrob(
    grid::gList(
      element_grob(calc_element("strip.background", theme_bw())),
      text_grob
    ),
    height = grid::grobHeight(text_grob), 
    width = grid::unit(1, "npc")
  )
}

# Plot digit function
scaleFUN <- function(x) sprintf("%.3f", x)

####Data####
# Setup working directory
WD = getwd()
setwd(WD)

respondents<-read.csv("respondents.csv")
cities<-read.csv("cities.csv")

####Figure 1####
# map
usmap<-map_data("usa")
statemap <- map_data("state")
cities$answer1<-as.factor(cities$answer)
cities$answer1<-ifelse(cities$answer==1,"Responded to Survey","No Response")

map.plot<-ggplot(usmap, aes(x=long, y=lat)) +
  geom_polygon(data = statemap, aes(x = long, y = lat, group = group), 
               colour = "black", fill = NA) +
  coord_map(xlim = c(-125, -65), ylim = c(24.5, 50)) +
  geom_point(data=cities, 
             aes(x=lng, y=lat, color=answer1),
             position=position_jitter(h=0.1, w=0.1),
             alpha = 0.5,size = 2) +
  guides(color = guide_legend(reverse = TRUE))+
  theme_article()

map.plot<-map.plot + theme(legend.position = c(0.2, 0.1),
                           legend.text = element_text(size = 12),
                           legend.title = element_blank(),
                           legend.background = element_blank(),
                           axis.title = element_blank(),
                           axis.text = element_blank(),
                           axis.ticks = element_blank())

map.plot
ggsave(file = "figure1.pdf", width = 8, height = 4, dpi = 666)
ggsave(file = "figure1.jpg", width = 8, height = 4, dpi = 666)

####Appendix A####
# descriptive stat summary
library(xtable) # function label will be masked
# city level
cities1<-subset(cities,cities$answer==1)
cityvars<-data.frame(cities1$population1000,
                     cities1$income1000,
                     cities1$leader_f_ratio,
                     cities1$labor_force_participation,
                     cities1$home1000,
                     cities1$unemployment_rate,
                     cities1$race_white,
                     cities1$race_black)

Mean=sapply(cityvars, mean,na.rm = TRUE)
SD=sapply(cityvars, sd,na.rm = TRUE)
Min=sapply(cityvars, min,na.rm = TRUE)
Max=sapply(cityvars, max,na.rm = TRUE)

city_des_stat<-matrix(NA,length(Mean),4)
colnames(city_des_stat)<-c("Mean","SD","Min","Max")
rownames(city_des_stat)<-c("Population","Median houshold income", 
                           "Female official ratio","Labor force participation",
                           "Home value","Unemployment rate",
                           "White percentage", "Black percentage")
city_des_stat[,1]<-Mean
city_des_stat[,2]<-SD
city_des_stat[,3]<-Min
city_des_stat[,4]<-Max
xtable(city_des_stat, digits=2)

# individual level
respondents$white<-ifelse(respondents$Race==1,1,0)
respondents$black<-ifelse(respondents$Race==2,1,0)
respondents$his<-ifelse(respondents$Race==3,1,0)
respondents$asian<-ifelse(respondents$Race==4,1,0)
respondents$other<-ifelse(respondents$Race==5,1,0)
respondents$female<-ifelse(respondents$Sex==2,1,0)
respondents$dem<-ifelse(respondents$Party==1,1,0)
respondents$rep<-ifelse(respondents$Party==2,1,0)
respondents$ideo<-ifelse(respondents$Ideology==6,NA,respondents$Ideology)
respondents$grad<-ifelse(respondents$Education>5,1,0)

vars<-data.frame(respondents$dem,
                 respondents$rep,
                 respondents$ideo,
                 respondents$Tenure,
                 respondents$white,
                 respondents$black,
                 respondents$his,
                 respondents$asian,
                 respondents$other,
                 respondents$female,
                 respondents$Age,
                 respondents$grad)
Mean1=sapply(vars, mean,na.rm = TRUE)
SD1=sapply(vars, sd,na.rm = TRUE)
Min1=sapply(vars, min,na.rm = TRUE)
Max1=sapply(vars, max,na.rm = TRUE)

des_stat<-matrix(NA,length(Mean1),4)
colnames(des_stat)<-c("Mean","SD","Min","Max")
rownames(des_stat)<-c("Democrats","Republicans","Ideology","Tenure",
                      "White","Black","Hispanic", "Asian","Other",
                      "Female","Age","Grad School")
des_stat[,1]<-Mean1
des_stat[,2]<-SD1
des_stat[,3]<-Min1
des_stat[,4]<-Max1
xtable(des_stat, digits=2)

####Appendix B####
#####Figure B.1####
# population
pop.plot<-cities %>%
  ggplot(aes(x=population1000, color=answer1)) +
  geom_density() +
  scale_x_continuous(trans = 'log10',
                     breaks = c(30,100,1000),
                     labels = c("30k","100k","1m"))+
  scale_y_continuous(labels=scaleFUN) +
  labs(x = "Population (Logarithmic Scale)",
       y = element_blank())+
  theme_bw()
pop.plot<-pop.plot+theme(legend.position = c(0.85, 0.90),
                         legend.title = element_blank(),
                         legend.background = element_blank())
pop.plot

# household income
inc.plot<-cities %>%
  ggplot(aes(x=income1000, color=answer1)) +
  geom_density() +
  scale_x_continuous(#trans = 'log10',
    breaks = c(50,100,200),
    labels = c("50k","100k","200k"))+
  scale_y_continuous(labels=scaleFUN) +
  labs(x = "Median Household Income",
       y = element_blank())+
  theme_bw()
inc.plot<-inc.plot+theme(legend.position = c(0.85, 0.90),
                         legend.title = element_blank(),
                         legend.background = element_blank())
inc.plot

# home value
home.plot<-cities %>%
  ggplot(aes(x=home1000, color=answer1)) +
  geom_density() +
  scale_x_continuous(trans = 'log10',
                     breaks = c(30,200,1000),
                     labels = c("30k","200k","1m"))+
  scale_y_continuous(labels=scaleFUN) +
  labs(x = "Home Value (Logarithmic Scale)",
       y = element_blank())+
  theme_bw()
home.plot<-home.plot+theme(legend.position = c(0.85, 0.90),
                           legend.title = element_blank(),
                           legend.background = element_blank())
home.plot

# labor
labor.plot<-cities %>%
  ggplot(aes(x=labor_force_participation, color=answer1)) +
  geom_density() +
  scale_x_continuous(#trans = 'log10',
    breaks = c(45,60,80),
    labels = c("45%","60%","80%"))+
  scale_y_continuous(labels=scaleFUN) +
  labs(x = "Labor Force Participation",
       y = element_blank())+
  theme_bw()
labor.plot<-labor.plot+theme(legend.position = c(0.85, 0.90),
                             legend.title = element_blank(),
                             legend.background = element_blank())
labor.plot

# unemployment
unemployment.plot<-cities %>%
  ggplot(aes(x=unemployment_rate, color=answer1)) +
  geom_density() +
  scale_x_continuous(#trans = 'log10',
    breaks = c(5,10,15),
    labels = c("5%","10%","15%"))+
  scale_y_continuous(labels=scaleFUN) +
  labs(x = "Unemployment Rate",
       y = element_blank())+
  theme_bw()
unemployment.plot<-unemployment.plot+theme(legend.position = c(0.85, 0.90),
                                           legend.title = element_blank(),
                                           legend.background = element_blank())
unemployment.plot

# female leader
f.plot<-cities %>%
  ggplot(aes(x=leader_f_ratio, color=answer1)) +
  geom_density() +
  scale_x_continuous(#trans = 'log10',
    breaks = c(0,40,100),
    labels = c("0%","40%","100%"))+
  scale_y_continuous(labels=scaleFUN) +
  labs(x = "Female Official Ratio",
       y = element_blank())+
  theme_bw()
f.plot<-f.plot+theme(legend.position = c(0.85, 0.90),
                     legend.title = element_blank(),
                     legend.background = element_blank())
f.plot

# white percent
w.plot<-cities %>%
  ggplot(aes(x=race_white, color=answer1)) +
  geom_density() +
  scale_x_continuous(#trans = 'log10',
    breaks = c(0,50,90),
    labels = c("0%","50%","90%"))+
  scale_y_continuous(labels=scaleFUN) +
  labs(x = "White Population",
       y = element_blank())+
  theme_bw()
w.plot<-w.plot+theme(legend.position = c(0.15, 0.90),
                     legend.title = element_blank(),
                     legend.background = element_blank())
w.plot

# black percent
b.plot<-cities %>%
  ggplot(aes(x=race_black, color=answer1)) +
  geom_density() +
  scale_x_continuous(#trans = 'log10',
    breaks = c(0,50,90),
    labels = c("0%","50%","90%"))+
  scale_y_continuous(labels=scaleFUN) +
  labs(x = "Black Population",
       y = element_blank())+
  theme_bw()
b.plot<-b.plot+theme(legend.position = c(0.85, 0.90),
                     legend.title = element_blank(),
                     legend.background = element_blank())
b.plot

# combine plot
represent.plot<-grid_arrange_shared_legend(pop.plot,
                                           inc.plot,
                                           f.plot,
                                           labor.plot,
                                           home.plot,
                                           unemployment.plot,
                                           w.plot,
                                           b.plot,
                                           ncol = 2, nrow = 4)

represent.plot<-annotate_figure(represent.plot,
                                left = text_grob("Density (%)",rot = 90))

represent.plot<-plot_grid(represent.plot)
represent.plot
ggsave(file = "figureb1.eps", width = 8, height = 10, dpi = 666)
ggsave(file = "figureb1.jpg", width = 8, height = 10, dpi = 666)

#####Table B.1####
popt<-t.test(cities$population1000[cities$answer==1],
             cities$population1000[cities$answer==0])

fratiot<-t.test(cities$leader_f_ratio[cities$answer==1],
                cities$leader_f_ratio[cities$answer==0])

inct<-t.test(cities$income1000[cities$answer==1],
             cities$income1000[cities$answer==0])

whtt<-t.test(cities$race_white[cities$answer==1],
             cities$race_white[cities$answer==0])

blkt<-t.test(cities$race_black[cities$answer==1],
             cities$race_black[cities$answer==0])

unt<-t.test(cities$unemployment_rate[cities$answer==1],
            cities$unemployment_rate[cities$answer==0])

homet<-t.test(cities$home1000[cities$answer==1],
              cities$home1000[cities$answer==0])

labort<-t.test(cities$labor_force_participation[cities$answer==1],
               cities$labor_force_participation[cities$answer==0])

repvar<-matrix(NA,8,3)
colnames(repvar)<-c("Responded Cities","No Response Cities", "P-value")
rownames(repvar)<-c("Population","Median houshold income", 
                    "Female official ratio","Labor force participation",
                    "Home value","Unemployment rate",
                    "White percentage", "Black percentage")

repvar[1,]<-c(popt$estimate[1],popt$estimate[2],popt$p.value)
repvar[2,]<-c(inct$estimate[1],inct$estimate[2],inct$p.value)
repvar[3,]<-c(fratiot$estimate[1],fratiot$estimate[2],fratiot$p.value)
repvar[4,]<-c(labort$estimate[1],labort$estimate[2],labort$p.value)
repvar[5,]<-c(homet$estimate[1],homet$estimate[2],homet$p.value)
repvar[6,]<-c(unt$estimate[1],unt$estimate[2],unt$p.value)
repvar[7,]<-c(whtt$estimate[1],whtt$estimate[2],whtt$p.value)
repvar[8,]<-c(blkt$estimate[1],blkt$estimate[2],blkt$p.value)

xtable(repvar, digits=2)


